
*****************************************
** Prep file that we'll use to make the tablepath

clear
	set obs 1 
	gen X = .
	save $datapath/temp/thetas.dta, replace
	
*****************************************


use $datapath/clean/main_data.dta, clear 
	keep if year<2020

*****************************************************************		
** Create the sigma and delta variables from demand estimation

	reshape long Q N, i(popu_lsa year lno ) j(type) s
	
	gen double M = 50*popu_lsa
	gen double s = Q/(N*M) 
	
	drop if s==.
	keep if s>0 & s<1 
	
	gen double de = type=="e"
	egen double QQ = sum(Q), by(lno year)
	
	gen double ls = ln(s) - ln((M - QQ)/M)
 
	gen double lstop = ln(Q/QQ)
	egen double NN = sum(N), by(lno year)

	gen double lstop_inst = ln(N/NN)
	gen double lsbottom = ln(1/N)
	gen double lsoverall = ln(1/NN)
	

	ivreghdfe ls  (lstop =   lstop_inst ) lsbottom   i.year de, absorb(lno) 

	gen double stop = _b[lstop]
	gen double sbottom =_b[lsbottom]
	
	
	gen double del = ls - stop*lstop - sbottom*lsbottom   
*****************************************************************		
	

keep if year==2018 
	
	
*****************************************************************		
** prepare variables (derivatives) needed to calculate the thetas
	
	gen double Dg = (N*exp(del/(1-sbottom)))^((1-sbottom)/(1-stop))
	egen double sDg=sum(Dg), by(lno year)
	gen double shat = (1/N)*Dg*(sDg^(-stop))/(1+sDg^(1-stop))
	
	gen double N1p = N + .1*(1-de)
	gen double N1e = N + .1*de 
	
	gen double Dg_1p = (N1p*exp(del/(1-sbottom)))^((1-sbottom)/(1-stop))
	egen double sDg_1p=sum(Dg_1p), by(lno year)
	gen double shat_1p = (1/N1p)*Dg_1p*(sDg_1p^(-stop))/(1+sDg_1p^(1-stop))

	gen double Dg_1e = (N1e*exp(del/(1-sbottom)))^((1-sbottom)/(1-stop))
	egen double sDg_1e=sum(Dg_1e), by(lno year)
	gen double shat_1e = (1/N1e)*Dg_1e*(sDg_1e^(-stop))/(1+sDg_1e^(1-stop))

	
	* note: dqp is the derivs of the two q's wrt p 
	* the wrt variable comes before the last, so - say - dqpe is the change in qe wrt np 
	gen double dqp = M*(N1p*shat_1p - N*shat)/.1 
	gen double dqe = M*(N1e*shat_1e - N*shat)/.1 

*****************************************************************		
	


tempfile loopdata
save `loopdata'


foreach g in none holdings {

	use `loopdata', clear 

*****************************************************************		
** Set prices 
	** materials costs calculated as mean expenditure/mean holdings from table 1 
	** staffing costs from section31_staffing_costs.do

	gen mp = 0.91
	gen me = 0.61
	
	gen g = "`g'"
	gen sp = 0 if g=="none"
	replace sp = 6.66 if g=="holdings"
	
	gen se = 0 if g=="none"
	replace se = 1.17 if g=="holdings"
	
	di "`g'" 
	su sp se
	
	gen pricep = mp + sp
	gen pricee = me + se
	
	gen pp = pricep
	gen pe = pricee

*****************************************************************		
	
	
	
*****************************************************************		
** Calculate thetas based on sigmas (and relative prices)

	preserve 
	
		keep lno type dqp dqe pp pe 
		reshape wide dqp dqe, i(lno ) j(type) string
		rename  dqpp dqp_dnp  
		rename  dqee dqe_dne  
		rename  dqep dqp_dne  
		rename  dqpe dqe_dnp  
		
		gen double det = dqp_dnp*dqe_dne - dqp_dne*dqe_dnp 
		
		gen double thetap = (dqe_dne/det)*pp - (dqe_dnp/det)*pe 
		gen double thetae = -(dqp_dne/det)*pp + (dqp_dnp/det)*pe 
		
		keep lno theta*
		reshape long theta, i(lno) j(type) string 
		drop if theta==. 
		egen double nn=count(theta), by(lno)
		keep if nn==2 
		save theta.dta , replace 
	restore 

	merge 1:1 lno type using theta.dta  
	keep if _merge==3
*****************************************************************		
	


*****************************************************************		
** Turn it into a wide dataset 	
	gen price = pp
	replace price= pe if type=="e"
	
	gen m = mp
	replace m = me if type=="e"
	gen staff = sp
	replace staff = se if type=="e"
	

	keep N theta stop sbottom price m staff M del lno type Q 
	
	
	reshape wide N theta  price m staff  del Q, i(lno stop sbottom M) j(type) string 
	
	gen ratio = thetae/thetap 
	
	keep lno ratio 
	gen g = "`g'"
	append using $datapath/temp/thetas.dta
	drop if g==""
	save $datapath/temp/thetas.dta, replace 
*****************************************************************		

}

*****************************************************************		
** Summarize thetas 
	
use $datapath/temp/thetas.dta
	
	su ratio if g=="none", de 
	su ratio if g=="holdings", de 
	
	count if g=="none" & ratio>1 & ratio!=.
	count if g=="none" & ratio<=1 & ratio!=.
